; FFT program by Piotr Nowak
; 10/21/2009
; Based on FFT C code by Tom Roberts, found on http://www.jjj.de/fft/fftpage.html



; fr[] and fi[] are real and imaginary arrays for input and output
; size of each array is 2^m

; mem(1000) = m
; mem(1001) = fr[0]
; mem(1001+2^m) = fi[0]
; mem(1001+2*2^m) = sinewave
addi r16, r0,0
bez r0, start

FIX_MPY:	; Fixed-point multiplication subroutine
		; parameters passed in r12 and r13, returned in r12
		; r13, r14 and r15 will be used for temporaties
		; r6 contains return address
mul r14, r15, r12, r13	; r14 = lo, r15 = hi
sll r12, r15, 1	
srl r14, r14, 16
srl r14, r14, 14
addi r15, r0, 1
and r13, r15, r14
add r12, r12, r13
srl r14, r14, 1
add r12, r12, r14

jalr r0, r6	; return

start:
li r1, 3
li r1, 232	; 1000, base mem address

; Find n = 2^m
load r2, 0(r1)	; r2 = m
addi r3, r0, 1	; r3 = 1
add  r4, r2, r0	; r4 = temp = m
nloop:
  bez  r4, nloopdone
  sll  r3, r3, 1
  addi r4, r4, -1
  bez  r0, nloop;
nloopdone:	; now r3 = n

addi r4, r0, 0	; r4 = mr = 0
addi r5, r3, -1	; r5 = nn = n - 1

; do decimation-in-time

addi r2, r0, 1 ; r2 = m = 1
ditForLoop:
  sub r14, r5, r2	; r14 = nn - m
  blz r14, endDitForLoop
  add r6, r3, r0	; r6 = l = n
  
  doLoop1:
    sra r6, r6, 1	; l >>= 1
    add r14, r4, r6	; r14 = mr+l
    sub r14, r5, r14	; r14 = nn - (mr+l)
    blz r14, doLoop1	; while( mr+l > nn )

  addi r14, r6, -1
  and  r14, r14, r4
  add  r4, r14, r6	; mr = (mr & (l-1)) + l  

  sub r14, r4, r2	; r14 = mr - m
  blez r14, contditForLoop ; if(mr <= m) continue;
	

  addi r7, r1, 1	; r7 = fr
  add r8, r7, r3	; r8 = fi

  add r14, r7, r2	; r14 = &fr[m]	
  load r13, 0(r14)	; r13 = fr[m]
  add r12, r7, r4	; r12 = &fr[mr]
  load r11, 0(r12)	; r11 = fr[mr]
  store r11, 0(r14)	; fr[m] = fr[mr]
  store r13, 0(r12)	; fr[mr] = fr[m]

  add r14, r8, r2	; r14 = &fi[m]	
  load r13, 0(r14)	; r13 = ti = fi[m]
  add r12, r8, r4	; r12 = &fi[mr]
  load r11, 0(r12)	; r11 = fi[mr]
  store r11, 0(r14)
  store r13, 0(r12)

  contditForLoop:
  addi r2, r2, 1	; ++m
  bez r0, ditForLoop	; end ditForLoop
endDitForLoop:	
	
; recycle r4, r5

addi r6, r0, 1		; l = 1
load r4, 0(r1)
addi r4, r4, -1		; r4 = k = m - 1

whileMain:
  sub r14, r3, r6
  blez r14, endWhileMain ; while (l < n) {


  sll  r5, r6, 1	; r5 = istep = l << 1;
  
  add r2, r0, r0	; m = 0
  outerForLoop:		; for (m=0; m < l; ++m)
    sub r14, r6, r2
	blez r14, endOuterForLoop


    add  r9, r2, r0	; r9 = j = m
    add  r14, r4, r0	; r14 = temp = k
    n2loop:
      bez  r14, n2loopdone
      sll  r9, r9, 1
      addi r14, r14, -1
      bez  r0, n2loop;
    n2loopdone:	; r9 = j = m << k;
    
    ; load twiddle factors
    add r12, r8, r3	; r12 = sinewave = fi + n
    sra r14, r3, 2	; n/4
    add r14, r9, r14	; j + n/4
    add r14, r12, r14	; &sinewave[j+n/4]
    load r10, 0(r14)	; r10 = wr = sinewave[j+n/4]
    add r14, r12, r9	; &sinewave[j]
    load r11, 0(r14)	; r11 = wi = sinewave[j]
    sub r11, r0, r11	; wi = -wi
    sra r10, r10, 1	; wr >>= 1
    sra r11, r11, 1	; wi >>= 1

    ; free up some registers, forget r1
    store r4, 3(r0)	; store k
    store r5, 4(r0)	; store istep
    store r3, 2(r0)	; store n
    store r6, 5(r0)	; store l
      

    add r1, r2, r0	; r1 = i = m
    store r2, 1(r0)	; store m
    innerForLoop:
      load r3, 2(r0)	; load n
      sub r15, r3, r1	; n - i
	  blez r15, endOuterForLoop  ; for (i=m; i<n; i+=istep) {

      
      load r6, 5(r0)	; load l
      add r9, r1, r6;	; j = i + l

      add r12, r10, r0	; r12 = wr
      add r15, r7, r9	; &fr[j]
      load r13, 0(r15)	; r13 = fr[j]
      addi r6, r0, 1	
      jalr r6, r6	; call FIX_MPY
      add r2, r12, r0	; tr = FIX_MPY(wr,fr[j])
 
      add r12, r11, r0	; r12 = wi
      add r15, r8, r9	; &fi[j]
      load r13, 0(r15)	; r13 = fi[j]
      addi r6, r0, 1	
      jalr r6, r6	; call FIX_MPY
      sub r2, r2, r12	; tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]);

      add r12, r10, r0	; r12 = wr
      add r15, r8, r9	; &fi[j]
      load r13, 0(r15)	; r13 = fi[j]
      addi r6, r0, 1	
      jalr r6, r6	; call FIX_MPY
      add r3, r12, r0	; ti = FIX_MPY(wr,fi[j])
 
      add r12, r11, r0	; r12 = wi
      add r15, r7, r9	; &fr[j]
      load r13, 0(r15)	; r13 = fr[j]
      addi r6, r0, 1	
      jalr r6, r6	; call FIX_MPY
      add r3, r3, r12	; ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fi[j]);

      add r15, r7, r1	; &fr[i]
      load r4, 0(r15)	; qr = fr[i]
      add r15, r8, r1	; 
      load r5, 0(r15)	; qi = fi[i]
      sra r4, r4, 1	; qr >>= 1
      sra r5, r5, 1	; qi >>= 1

      sub r14, r4, r2
      add r15, r7, r9	      
      store r14, 0(r15)	; fr[j] = qr - tr

      sub r14, r5, r3
      add r15, r8, r9	      
      store r14, 0(r15)	; fi[j] = qi - ti


      add r14, r4, r2
      add r15, r7, r1	      
      store r14, 0(r15)	; fr[i] = qr + tr


      add r14, r5, r3
      add r15, r8, r1	      
      store r14, 0(r15)	; fi[i] = qi + ti


      load r5, 4(r0)	; load istep
      add r1, r1, r5	; i += istep
      bez r0, innerForLoop
    endInnerForLoop:
      
    load r4, 3(r0)	; load k
    load r5, 4(r0)	; load istep
    load r3, 2(r0)	; load n
    load r6, 5(r0)	; load l
    load r2, 1(r0)	; load m
    addi r2, r2, 1	; m++
    bez r0, outerForLoop
  endOuterForLoop:

  addi r4, r4, -1	; --k
  addi r6, r5, 0	; l = istep
  bez r0, whileMain

endWhileMain:

